

shootings <- readRDS("temp/geocoded_shootings_from_gva.rds") |> 
  filter((date >= "2020-05-03" & date <= "2021-05-02") |
           (date >= "2016-05-08" & date <= "2017-05-07"),
         vics >= 4) |> 
  mutate(post = (date >= "2016-11-08" & year(date) %in% c(2016, 2017)) |
           date >= "2020-11-03",
         year = 2 * floor(year(date) / 2),
         d2 = ifelse(year == "2020", as.integer(as.Date(date) - as.Date("2020-11-03")),
                as.integer(as.Date(date) - as.Date("2016-11-08"))))


gangs <- rdrobust(y = shootings$gang, x = shootings$d2, p = 1, c = 0)
dvs <- rdrobust(y = shootings$dv, x = shootings$d2, p = 1, c = 0)
home <- rdrobust(y = shootings$home, x = shootings$d2, p = 1, c = 0)
deaths <- rdrobust(y = shootings$deaths, x = shootings$d2, p = 1, c = 0)
vics <- rdrobust(y = shootings$vics, x = shootings$d2, p = 1, c = 0)


outs <- rbindlist(lapply(c("gang", "dv", "home", "deaths", "vics"), function(v){
  
  l <- rdrobust(y = shootings[[v]], x = shootings$d2, p = 1, c = 0)
  
  tibble(var = v,
         coef = l$coef,
         u = l[["ci"]][,2],
         l = l[["ci"]][,1])
})) |> 
  mutate(across(c(coef, u, l), ~ . * -1))

outs$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(outs)/3)

outs <- filter(outs, estimate == "Robust") |> 
  mutate(var = case_when(var == "gang" ~ "Gang Related",
                         var == "dv" ~ "Domestic Violence",
                         var == "home" ~ "Home Invasion",
                         var == "deaths" ~ "Number Killed",
                         var == "vics" ~ "Number Injured"))
## figure s18
ggplot(outs, aes(x = var, y = coef, ymin = l, ymax = u)) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 10,
           plot.margin = margin()) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  # scale_x_continuous(breaks = seq(0, 10, 2)) +
  # scale_y_continuous(limits = c(-0.1, 0.2)) +
  labs(y = "Local average treatment effect", x = "")

for(i in outpaths){
  ggsave(paste0(i,"shooting_characteristics.png"), height = 4, width = 5, units = "in")
}

dens <- rddensity(shootings$d2, c = 0, p = 1, h = 51.5)
summary(dens)
rdplotdensity(dens, shootings$d2)
